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ABSTRACT 

The partial differential equation for heat diffusion is numerically 
integrated by the Runge-Kutta-Gili method. A solution is obtained for 
the diurnal temperature variation with a bounded coefficient of eddy 
diffusivity which varies periodically with time and nonlinearly 
with height. The surface wave is represented by the sum of a 
diurnal and a semidiurnal harmonic wave. The results may be interpreted 
to apply over a fairly broad range of diffusivity values and height. 
With avpropriate choices of the various parameters, reasonably good 
agreement is obtained between theoretical and observational values of 
amplitude reduction and phase lag. 

The author wishes to express his appreciation for the assistance 
and encouragement of Professors G. J. Haltiner and E. J. Stewart of 


the U. S. Naval Postgraduate School in this investigation. 
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ei20rt nas peer directch tows. %e diurmot termeriturc wove, ondwia 
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por, tcular location. 

In many of the past studies, « has been asuyuwnel to vary, (a) 
dineaciy with heisht, (tb) as a vower of height, and (c) as a bounic: 
exponentiar Function of neiaht. In each of these cases, the essumed 
Meect ional Porm OF A was agsumed tc hold tirousnout the etruosohere. 
Recertly, nowever, de Vries (1) has studied: 


the case of a: atmosynnere consistin: of several levers, 
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eacn of tnese layers, ni in particular, the turbulert boundary 
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analytical solutions for several particular cases; (a) K = constant, 
(bd) K =Xtese,)”, and {c) K = ali - b exp(-cz)] . Case (a) yields a 
solution in terms of a rether complicated exponential function, (b), 
a solution in terms of Bessel functions of the first kind, and (c), 
a@ sOlution in terms of a hypergeometric function. However, no 
numerical values were given, and the results are not suitable for 
practical vurposes. 

It is the puroose of this investization to present a numerical 
solution of Eqn. (1), with appropriate boundary conditions, for a 


varticular case of an atmosphere composed of an arbitrary number of 


layers. 





2% The Heat Diffusivity Coefficient. 

From physical considerations, K, the heat diffusivity, may be 
expected to increase in the lower layers, but should not continue to 
increase indefinitely with height. Hence, a bounded diffusivity 
coefficient, which increases with height from the surface to some fixed 
level and then decreases with height to some residual value at higher 
levels, would seem more suitable than either the case of K increasing 
linearly with height or of K increasing as a simple power of height. 
Upon examination of Fig. (1), which is based on observations taken 
at the micrometeorological tower at the University of Washington (2), 
(3) and Mildner's pilot balloon observations taken near Leipzig (5), 
it is seen that no single functional form of K will suffice. The 
values of K from the University of Washington data were obtained 
by methods based on energy continuity at the earth's surface, 
while those from the Leipzig data were actually values of eddy 
viscosity. These values of eddy viscosity were used as an 
approximation to eddy diffusivity, since they showed the desired 
variation with height, and since values of eddy diffusivity at 
sufficiently great heights are difficult to obtain. Upon 
dividing the atmosphere into four layers, and then fitting functions 


of height to each layer, the following form of Kk resulted: [see Fig. 


(2) 
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The quantities By Bg yB 5184s Oe 8, and h are constants, appropriate 
values of which will be assigned later. The discontinuity which 
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occurs at z = 228 meters is resolved by defining FR at this level as zero. 
Haltiner (4) suggests that K should have a diurnal variation, and that 
it appeers reasonable to expect that K woulec more or less follow 


the surface temocrature wave, reachings maximum and minimum values 
near the time when the temperature is maximm and minimm, 


respectively. 
The eddy diffusivity, K, can now be expressed as a function of height times 
a function of time, i. e., K = f(z)e(t), where g(t) is defined: 

a(t) =1 + b(sinWt + c¢ sin ewt) (3) 
Here b and c are constants, and W «= 277/86,400 Sec. 

Appropriate boundry conditions associated with Eqn. (1) will now be 

chosen. These are frequently taken to be: 

@ =0 for z =O, and all t; (4) 

&Y = O_ (t) for 2 2 0; (5) 
where O (+) is defined by: 

O(t) = dlsinwt + ¢ sin ewt) (6)- 


Here d is a constant and c and@are the same as in Eqn. (3). 








oP Transformation of Coordinates. 

In order to piace the basic equations into a form which will give 
a broader interpretation of the results, and to scale these equations 
for computational purposes, let: 

2 =100qgo, and t = 36007. (7) 

Here q is a constant, while © and & are the new variables. When t is 
in seconds, the units of “ are hours, and when z is in centimeters, the 
units of co are in "q-meters"; 1. e., in lemeter units when q = 1, 2- 
meter units when q = 2, etc. With this transformation, (7), Faqs. (1) 


and (2) become: 


J ee Vee L 


(3) 


| 2, ome 
Ki = Hel 4(¢) : q(t) Q'S +a; 


Gy Expl -% le -A 


(9) 


Ae 
Here WY must now be taken as 2/!//2h, and; 
a! sa, (100q)*, a't = a,(100q), 
at = a, (100q), h' = h/100q. 


The boundary conditions, Eqs. (4) and (5), remain identical in form. 
Both observation and theory indicate that the most pronounced 
variations of temperature occur near the surface of the earth. 
Tnerefore, 2 small grid distance is desirable at low levels, while 
at higher levels, a small grid distance is not necessarily needed. 
This suggests that the following transformation of the vertical 
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coordinate may be useful: 

s=ino. (10) 
With this transformatior, the lower boundary condition, Eqn. (5), may 
apoly at a heiztht of a meters corrosponding to © = 1, however, to reduce 
surface effects, it is convenient to select © = 2. Actually the vertical 
coordinate system is arbitrary, and can be chosen to have its base at any 
level. Py utilizing Man. (10), Tas. (8) and (9) now become: 


IG 25h a, 
o ee jy) te 
zB / ies : see: Sa +(Gh - K)F< | 





J (11) 


Q, ‘e° 
K = Ks) g(t) = GA(T) \ ae +a; 
hy x P| ~4% !(E>-4 | 
6 (12) 
Here avy is taken as On. and the other constants, a) ; aD) Bo, Oa, 


and ag, must be scaled oy a factor of 10 °. 





4. Finite Difference Equations. 

In order to obtain a numerical solution, the derivatives of @ 
with respect to s and @ in Eqn. (11) must be replaced by appropriate 
finite difference forms. The problem is then reduced to that of 
sOlving a system of linear algebraic equations in the values of O 
over a grid of voints covering the desired range. Expressing 3-2. 
and 9% in Fan. (11) as finite difference ratios, we obtain: 

d @. = 0. 7 a 04, -2& ag 
ae a4 Cc (M,) ( ae 


Pie 4) S23) 


L242 5% ---- mw. 
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Here the subscript 1 designates the i'th level of the vertical grid 
at height 12, where P is the vertical distance between grid points in 
units of s, that is, a pure number , Since s is dimensioniess. 

The partial differential equation, Eqn. (1), has now been reduced 
to a system of ordinary differential equations which will be 


integrated numerically by the Runge-Kutta-Gill method. 





oe The Constants. 
= ; i a, 2 
Lext, appropriate values off, on, anc AC, as well as the constants 


will now be chcser. ince the vertical cooréinate is 


Ap, Az, ete., 
s- inc, Ye S, a constanm®. With @ choice of & = in 2, thew 
equation of the systen defined by Eqn. (12) apnlies at the height 

e cd-Meters. Therefore, tne wuecessive valves Of 1, besinning with 

the lower boundary condition, apply to the levels, 2, 4, 8, 216, 32,5k, 
etc, q-teters. To meduce thewinve tration time, yet still ive «= tality 
eetailel nicture Of the vertical structure;eene Upoer boundary ccnditior 
"an. (4, 1s chosen to anply at the height of 1024 g-meters. Yor this 
choice, Ge =%; and the system, Eqn (13), to be solved consists of 
eight sivulteneous ordinary differential equations, together with 

riqn . (6) for ey A numerical solution for this system was then 
obteined by the Runge-“utta-%i1ll method on a National Cash Register 

LOCA electronic comouter. ‘The choice of an appropriate time interval 
Yepended intimately upon the size of the eddy diffusivity. As the 
Patver increased, At. had to be decrease? in order to maintain 
computational stability and reascnable accuracy. 


The following numerical values of the other constants were 


enosen as bein revresentative: 


ay = 9.% (sec os aS = 2050 (cm/sec), 

Be = =92,900 (en“/sec ), 2 Wi O00 (er*/sec), 

a} = 9.9085, a, = 58,000 (em*/sec), 

eee yo. ooe= C.K 

GeO. O54 Ce Nh’ = 220. (14) 


The cjoice of a! 


Peete earn G7, giver f 525-1010 increase 


of K with heignt from tue surface to 223 a-meters, and a 14-fold 
jecrease of K with heigtt above 22d q-metercs. ‘Also, the cnoice of 
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b and c gives rise to a 21-fold diurnal variation of K with time. The 


value of d was chosen because it keeps the magnitude of 9, <0.1. 

With the value of ay= 10° , the integration interval was 
prohibitively short for the computer used; therefore a smaller 
value of a7 = G-Poas 10* was used. This permittec a time interval 
of 1/30 hour. 


The value of ay = 6.25 x 104 gives the following range of the 


coefficient of eddy diffusivity, Eqn. (12), in units of cm*/sec: 


at surface, Kp 0-2 Sorc. K pip ae a 
at 20 q-m, Kmax = 0.15 3 10-ce: Ky fee f 14 A= 
at 228 q-n, Kmay = 0-492 x 1Ocqc, Seyets 2a alee 
at 460 q-n, K may = 0.35 x 10°q*, nope Ste ae 


For q £1, the values of K are somewhat small for typical etmospheric 
conditions, esvecialily in the lowest layers. This is reflected in 
the results by a large amplitude reduction and vhase lag with 
increasing height. However, such values of K are not unrealistic 


for winter. 
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The results of the interration are given »y Pi. (3) and Table ee) 
Fig. (%) being the graphical internretation of Table (1). Tr Pir. (5), 
Syplitude reduction snd phase tae, in hours, are pictted apainst Aegis 
in q-meters. As enn be seen from Fir. (3), the amolitude reduction and 
Dnase lar are quite larsc, especially wnen compared with ovservei values. 
(See Fig. (4) which is besec on observations taken at the micrometeorological 
tower at tne University of Washinton). The large armlitude reduction 
ans phase las from the theoretical results, Fig. (3), are due largely 
to the reduction of a from 10° to “~.25 x 10%. however, by vrover 
@ecection Of tdervalic OF Gy fai Serecmemranecyees Ticoret tea. and 
observed values can ve obtained. Fig. (4%) also shows the amplitude 
reduction and nhase lag for the mininuw: at three levels. From these 
tiree levels, it can be seen that the anmlitude reduction is less 
moi the maximum tman for the minimum, che tatvcer corresvonding to the 
lower nocturnal 7fusivity. This difference is anpreciatie in that at 
the eitht q-meter level, the mxinmum ts about 43°%/o of the surface value 
ohm © white the minimum is only 19°/o of the surface minimum. At lower 
Hevels, however, Unis difference is auch samiller. The phase law increases 
with heisht for both the meximum and the minimum, however, the laz 
of tre maximum is greater than that of the minimum at any varticular 
jJevel. “his does not seen consistent with greater amplitude reduction of 
minina compare? to maxima. Possibly, these cifferences would vary 
mo, SucceeJing waxi-1a and minima, as preater time elapses from the 
madi ial concitions. 

By a proper choice of q, some comparisons vetween observed and 
theoretical results may be mde. Table (2) shows amolitude reduction 
anc ohase lag from cbservations taken at IToafield, alonp with 
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theoretical values for q = 10 and q = 20. For voth values of q, there 

is fair agreement with observed values, however, both the theoretical 
amplitude end phase lag increase with height much more rapidly than do 
the observed values. Table (3) comoares theory with observed data at 

the Eiffel Tower. Here q is taken to be 60 and 33, the value of 60 
comparing to the summer data and the 33, to the winter data. Table (4) 
shows some observations taken at the University of Washington. By the 
choice of q = 8, the amplitude compares fairly well; but where the observed 
phase lags are the order of minutes, those from theory are the order 

of hours. In Tabdie (5), ais ei from Porton is comoared to theory. For 
the choice of q = 2.5, the comparison is fair, but once again, shows 

the theoretical phase lag anc amplitude reduction much greater than 

the observed. Table (6) commarcs the amplitude reduction in summer 

and winter at Teafield with theory. Here q = 5 was chosen for winter, 
and q = 12, for summer. For the winter case the theoretical and 
Observed values compare very well, but for summer, the theoretical 


amplitude reduction is simmificantly greater than the observed value. 


“Fe 


1 





even Tg east > Ld : 
Bepiicude regicticn (2. F. 


ent AT) «% By ~ p> A 
7 ee LOG q eee Os ous 13 caahee ha Qiise a 





rd 
adnan eens TV et ee c 2 
Peewee CUO ree Of Fro ieee ieee i tee ve sors CIeOreyrer! 


V44rms! eae : — ; S 
Poe UC £OY oa = oe anc 0 mee 


I Oe EELS EG, © TLE AN EPO 





a eee eed 


lira / Wet Rice 83 











PabtC ya. 
Amplitude reduction and phase lag for tne case of gq = & versus 


observations taken at the University of Washington. 
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Aymmlitude reduction and phase lag for the case of q = @.5 versus 


Ooservations taken at Porton, “ingland. 
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Amlitude reduction from tue leafiels data versus theoretical 


Beaiuecs for gq = 5 ard 1 = 12. 
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i= Conclusions. 

A numerical solution has been presented for the diurnal temperature 
variation with a coefficient of eddy diffusivity which is a function 
of height and time. By suitable selection of several parameters, 
reasonably good agreement has been obtained between theory and 
observation. 

Improvement in the results may possibly be achieved by improving 
the functional form of the diffusivity coeffictent, especially in the 
lowest layers. It would also be desirable to include such parameters as 
surface roughness, stability, wind velocity, etc. With a diffusion 
eserficient im Pri of these parameters, as well as height and time, 

@ sOlution for the diurnal temperature wave would have a mich broader 


application. 
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